ANALYSIS OF A DARCY-CAHN-HILLIARD DIFFUSE INTERFACE 
MODEL FOR THE HELE-SHAW FLOW AND ITS FULLY DISCRETE 
FINITE ELEMENT APPROXIMATION 

XIAOBING FENG* AND STEVEN WISE"!" 

Abstract. In this paper we present PDE and finite element analyses for a system of partial 
differential equations (PDEs) consisting of the Darcy equation and the Cahn-Hilliard equation, which 
arises as a diffuse interface model for the two phase Hele-Shaw flow. In the model the two sets of 
equations are coupled through an extra phase induced force term in the Darcy equations and a fluid 
induced transport term in the Cahn-Hilliard equation. We propose a fully discrete implicit finite 
element method for approximating the PDE system, which consists of the implicit Euler method 
combined with a convex splitting energy strategy for the temporal discretization, the standard finite 
element discretization for the pressure and a split (or mixed) finite element discretization for the 
fourth order Cahn-Hilliard equation. It is shown that the proposed numerical method satisfies a 
mass conservation law in addition to a discrete energy law that mimics the basic energy law for the 
Darcy-Cahn-Hilliard phase field model and holds uniformly in the phase field parameter e. With help 
of the discrete energy law, we first prove that the fully discrete finite method is unconditionally energy 
stable and uniquely solvable at each time step. We then show that, using the compactness method, 
the finite element solution has an accumulation point that is a weak solution of the PDE system. 
As a result, the convergence result also provides a constructive proof of the existence of global-in- 
time weak solutions to the Darcy-Cahn-Hilliard phase field model in both two and three dimensions. 
Finally, we propose a nonlinear multigrid iterative algorithm to solve the finite element equations 
at each time step. Numerical experiments based on the overall solution method of combining the 
proposed finite element discretization and a nonlinear multigrid solver are presented to validate the 
theoretical results and to show the effectiveness of the proposed fully discrete finite element method 
for approximating the Darcy-Cahn-Hilliard phase field model. 

Key words. Two phase Hele-Shaw fiow, diffuse interface model, Darcy law, Cahn-Hilliard 
equation, energy splitting, finite element method, nonlinear multigrid. 
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1. Introduction. Hclc-Shaw flow refers to the motion of (one or more) viscous 
fluids between two flat paraUel plates separated by an infinitesimaUy smaU gap. Such 
a physical setup is often called a Hele-Shaw cell and was originally designed by Hele- 
Shaw to study two dimensional potential flows |17j . Various fluid mechanics problems 
can be approximated by Hele-Shaw flows and thus the research of those flows is 
of great practical importance. In addition, the relative simplicity of the governing 
equations of these flows makes Hele-Shaw flows ideal test cases in which rigorous 
mathematical theory and efficient numerical methods can be developed for studying 
interfacial dynamics — such as the formation of singularities and topological changes 
— in immiscible fluids (cf. [20l[2TJ|24] and the references therein). 

The governing equation of Hele-Shaw flow is identical to that of the inviscid 
potential flow and to the flow of fluids through a porous medium, because the gap- 
averaged velocity of the fluid is given by Darcy's law. Speciflcally, the two phase 
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Hele-Shaw flow takes the form (cf. [20] and the references therein): 

(1.1) u=-^(Vp-pg), divu = innrXTt, 

(1.2) [p] = -fK, [u • n] = on Ft, 

with a given set of initial and boundary conditions. Here fix = ^ ^ {0,T), where 
17 C is a bounded domain. Tt denotes the interface between the fluids at the 
time t with the normal n. u is the fluid velocity and p stands for the pressure of 
the fluids. The symbol \p] stands for the jump of p across the interface Ft- r] is 
the viscosity, which may have different (positive constant) values on both sides of 
Ft. g is the gravitational force per unit mass; and p is the mass density of the fluid, 
which again can take different (positive constant) values on both sides of the interface. 
Equation ( 1.1 )a is Darcy's law and ( |l.l[ )b implies that the fluids are incompressible. 



Equations (1.2)a and (1.2 lb are the boundary conditions at the fluid-fluid interface. 



which represent the mathematical descriptions of the balance of the surface tension 



forces and the balance of mass, respectively. Equation (1.2)a is called the Laplace- 
Young condition, where 7 is the dimensionless surface tension coefficient and k is the 
(mean) curvature of the interface Ft. Notice that the tangential component of the 
velocity u may experience a jump across the interface |20j . 

Computationally, the above moving interface problem is difficult to solve directly 
due to the existence of the surface tension on the interface. In addition, during the evo- 
lution the fluid interface may experience topological changes such as self-intersection, 
pinch-off, splitting, and fattening. When that happens, the classical solution of the 
moving interface problem ceases to exist. In such cases it is very delicate matter to 
develop a proper notion of generalized solutions, and it becomes even more challenging 
to compute the generalized solutions when they can be defined. 

To overcome the difhculties, an alternative approach for solving moving interface 
problems is the diffuse interface theory, which was originally developed as methodol- 
ogy for modeling and approximating solid-liquid phase transitions in which the effects 
of surface tension and non-equilibrium thermodynamic behavior may be important 
at the interface. In the theory, the interface is represented as a thin layer of finite 
thickness, as opposed to a sharp interface. Such an idea dates to Poisson, Gibbs, 
Rayleigh, van der Waals, and Korteweg (see [231 HO] and the references therein). The 
approach then uses an auxiliary function (called the phase field function) to indi- 
cate the "phase". The phase field function, denoted by (p below, assumes distinct 
values in the bulk phases away from the interfacial region, through which the phase 
function varies smoothly. The interface itself can be associated with an intermedi- 
ate contour/level set of the phase function (cf. [H |31 [51 [121 [53] and the references 
therein) . Generally speaking, the diffuse interface models are expected to converge to 
some corresponding sharp interface models as the width of the interfacial layer tends 
to zero. 

The diffuse interface model for Hele-Shaw flows to be studied in this paper is 
given as follows: 

(1.3) u = — Vp — 7(/3V^ in r^T, 

(1.4) div u = in ftr, 

(1.5) (ft + u ■ \7(p — eAp. = in r^T, 

(1.6) n=-eAip+^f{^) in flT, 

£ 
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where f{ip) — P'iv) ^nd F{ip) — \{(p'^ — 1)'^ is the so-called double-well (potential) 
energy density, and < e < 1 is a fixed constant. To close the system, we impose the 
following initial and boundary conditions 



(1.7) 
(1.8) 



dp 
dn 



9/1 dip 
dn dn 







da 



T •- 



dn X (o,r], 



n. 



Note that we have suppressed the superscript e in {u'^ jip'^) for the sake of nota- 
tional simplicity. Although $7 is a two-dimensional domain in the original Hele-Shaw 
problem, in this paper we consider C M'' {d — 2, 3) because the three-dimensional 
problem also has a mathematical interest and arises from biological applications |28j . 
Here the vector u{x,t) G M.'^ and the scalar p{x,t) E M denote the velocity and 
the pressure of the fluid mixture at the space-time point {x,t), respectively. The 
variables tp{x,t), fJ-{x,t) € K are the phase field function and the chemical poten- 
tial, respectively, ip assumes distinct values — namely, ±1 based on our choice of 
F(ip) — in the bulk phases away from a thin layer of width 0(e). This thin layer 
is called the diffuse interface region. It is natural to define the zero level curve of 
ip, T^{t) — e M'^ I ip{x,t) = O}, as the d — 1 dimensional interface. Eq. (1.3 1 with 
7 = is the Darcy equation [4] . (|1.5[) and ( 1.6 ) without the convection term u • V^s is 



the Cahn-Hilliard equation [8} 111 ! 123 ] . Note that if 7 = 0, the velocity vanishes, and 
the Cahn-Hilliard equation results. 

The system ( |1.3| -( |L81 ) is a special case of the BHSCH (Boussinesq-Hele-Shaw- 
Cahn-Hilliard) model proposed by Lee, Lowengrub, and Goodman in j20j . They 
showed, using formal asymptotics, that solutions of the BHSCH system converge to 



those of the Hele-Shaw model (1.1H(1.2) as their interfacial parameter converges to 



zero. We note that the pressure p in (1.3 1 has a different scaling from that in the 



as the 



BHSCH model in [20]. To obtain a similarly scaled pressure, one can simply introduce 
a redefined pressure in our model via p — p + "fffJ-- We shall refer ( 1.3 )-( 1.8 
DCH (Darcy-Cahn-Hilliard) model/system herein. 
Define the Cahn-Hilliard energy functional 



(1.9) 



e.^ ,2 1 



.\Vip\' + -Fiip) 



dx. 



Like many diffuse interface models (cf. [31 [T^ [231 UHl 123])) the DCH system is also 
a dissipative system as it satisfies the following energy dissipation law (see Sec. [2] for 
the details): 



(1.10) 



dt 



IVa^II 



1 

7 



= 0. 



As expected, the above energy law plays a vital role in the analysis of the DCH system 
and in the design and analysis of numerical methods for the system (see Sees. |2jj4] for 
the details). 

This paper consists of four additional sections. Section [2] is devoted to the PDE 
analysis of the initial-boundary value problem ( 1.3 1-(|1.8|). Weak solutions are de- 



fined and the uniqueness and regularities of weak solutions are established. Section |3] 
contains the formulation of our fully discrete implicit finite element method for prob- 



lem (1.3)-(1.8). It is shown that the proposed numerical method satisfies a mass 



conservation law in addition to a discrete energy law that mimics the basic energy 
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law for the Darcy-Cahn-Hilliard phase field model and holds uniformly in the phase 
field parameter e. With help of the discrete energy law, it also proved that the fully 
discrete finite method is unconditionally energy stable and uniquely solvable at each 
time step. Section |4] presents a convergence analysis for the proposed fully discrete 
finite element method. Using the compactness method it is shown that the finite 



element solution has an accumulation point that is a weak solution of problem (1.3) 



(1.8). As a byproduct, this convergence result also provides a constructive proof of 



the existence global-in-time weak solutions to the PDE system (1.3)-(1.8| in both 
two and three dimensions. Finally, in Sec. [5] we provide some results of numerical 
experiments validating our theoretical results and showing the effectiveness of the 
proposed fully discrete finite element method. To solve the nonlinear finite element 
equations at each time step, we propose a nonlinear multigrid iterative method to do 
the job. The details of the nonlinear multigrid solver and some other algorithmic and 
implementation issues are described in Appendix [A] 



2. PDE analysis. The standard space notations are used in this paper, we 
refer to [U [9] for their exact definitions. In particular, B* denotes the dual space 
of a Banach space B, and B denotes the vector Banach space B"^, where d is the 
dimension space. Here we shall assume d = 2 or 3. The symbol (•, •) is used to denote 
the standard inner product, (•,•) stands for the dual product between H^{rt) 

and {H^{n))* . ig(57) denotes the subspace of L^(f2) whose functions have zero mean. 
Throughout the paper, unless stated otherwise, c and C will be used to denote generic 
positive constants which are independent of p, /x, ip, u, and e. If, for example, there is 
a dependence on e, we shall explicitly write C — C{e). As indicated earlier, we shall 
assume that < e < 1. 

In the next section we shall construct a finite element method which directly 
approximates variables p, /i, and ip, but not u, which will be computed as an auxiliary 
variable as needed. Specifically, we shall approximate the pressure equation by a 
standard finite element method and the phase equation by a mixed finite element 
method. We remark that it is also a viable strategy that approximates both the 
pressure equation and the phase equation by mixed finite clement methods, which we 
shall study in a separate work. 

To get the governing equations without using u, substituting the expression of u 



in (1.3) into (1.4) and (1.5) we get 



(2.1) div (Vp 7</3V/i) = in r2r, 

(2.2) (pt-eA^-div ((^[Vp-f 7(/3V^]) =0 in f^T. 

Then the PDE system to be studied and approximated in this paper consists of 



initial conditions (1.7)- (1.8) 



equations (2.1), (2.2), and (1.6), which are complimented with the boundary and 



Motivated by the energy law ( 1.10 ), we define the following weak formulation and 



solutions to the initial-boundary value problem. 

Definition 2.1. Let (p^ e H^{VL). A triple {p,^,(p) is called a weak solution of 
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problem (2.1), (2.2), and (1.6)~(1.8) if it satisfies 



(2.3) peLl {io,Ty,H\n)nLl{n)), 

(2.4) ^eL2((o,r);i/i(r!)), 

(2.5) Vp + 7^VMGi'((0,r);L2(O)), 

(2.6) ipeL°" {iO,T);H\n)) n {{0,T); (n)) 

(2.7) ipteLl{{0,T);{H\n)r), 

and there hold for almost all t € (0, T) 



(2.8) 
(2.9) 

(2.10) 



(Vp + -fipVfi, Vq) = 
((^[Vp + 7(pV^] , Vi/) = 



e(V(p,V^) --(/(¥>), V^) = 

£ 



wii/i f/ie initial condition ip{0) — (^g. 

Remark 2.1. The reason for not breaking the sum Vp + "fipS/^ is that it has 
better regularity /integrability than each of its two terms. By the Aubin-Lions lemma 
(cf. !l25f ). the regularity on ip ensures that if e C"'([0,T];L2(f])). Hence, the initial 
condition (p(0) = (ySg makes sense. 

Remark 2.2. The regularities imposed on the solution (p,fi,ip) in the definition 
are not the "minimum" required to make all terms in (2.8)-(|2^T0) be well defined. 
These regularities are imposed because they are suggested by the energy law (1.10). 



Moreover, the product space of the spaces used in the definition for {p, fi, ip) is indeed 
the energy space associated with the DCH system. 

As we mentioned in Sec. [T] a key feature of the Darcy-Cahn-Hilliard (DCH) 
system is that it is a dissipative system in the sense that it satisfies an energy law, 



namely, (1.10). Below we demonstrate that this is indeed the case for weak solutions 



of the DCH system. 

Lemma 2.2. Suppose that (^g e H^{n) and Vl C M.^_{d = 2,3) is a Lipschitz 
domain, let (p,ii,ip) be a weak solution defined by (2.8)-( 2".10[ ). In addition, suppose 
that the initial value ip^ satisfies J,, ((/3g) < Cq for some e-independent constant Cq, 
i.e., the initial energy is uniformly bounded in e. Set u := — (Vp + jipVn). Then, for 
almost all t G (0,T), 



(2.11) 



iy9(a;, t) dx 



dx, 
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and for all t G (0,T) and some e-independent constant C — C(i?(0)) > 



(2.12) 

(2.13) 

(2.14) 

(2.15) 

(2.16) 

(2.17) 
(2.18) 



E{t) + / e \\'^^i{s)\\l. + - \\n{s)\\l, ds = £;(0) < oo, 
Jo L 7 -I 

max ||</'(s)||iri < — , 

Ms)\\l, ds < ^ , 



-e V(v(s))||r2 ds < 



(T+1)C 



\\ft{s) + u(s) • V¥3(s)||(^i). ds < Cs, 
C 



ds < 



max |||<y9(s)| - 1||^2 < Ce, 

[)<s<t 



where E(t) :— (}p(i)) and ■ ) is defined in (1.9). 

Proof. (2.11) follows trivially from setting = 1 in (2.9). To prove (2.12), we 
consider two cases separately. First, suppose that ipt G L'^ ((0 , T); H^{^)). Then 
setting g = 2 in (2.8), v — fj, in (2.9) and ^ — —ipt in ( 2.10| ) (note that by the 
assumption, —ipt is a valid test function), and adding the resulting equations we get 



d r 
dt 



e II Vy^lli. + -(F(^), 1) + e WV^Wi, + - ||Vp + -f^V^i\\i2 = 0. 



Integrating over the interval {0,t) yields (2.12). 

For the general case (ft G ((0, T); (if^(ri))*), and we note that V = — is not 
a valid test function in (2.10). However, this technical difficulty can be overcome by 
using a Steklov average technique. For t G (0, T), let 5 > be a small number. Define 
the Steklov average (p^ of ip by (cf. [HI Ch. 2]) 



^\-,t) :=4(^)(.,i) 



t+5 



ip{-,s)ds vte(o,r). 



Trivially, for small enough S, 



^f{-,t) {v'{;t))^ 



ipi-,t + d)-^{-,t) 



Hence, (Pt{ - ,t) € i7^(f2) for almost every t £ (0, T- (5). It is well known that (cf. [HI 
Ch. 2]) 



(2.19) 



Note that the derivative on the left-hand side of the above identity is understood as 
a distributional derivative, while the derivative on the right-hand side is understood 
in the classical sense. 
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Now applying to both sides of ( |2.8[ )-(2.10) after replacing i by s yields 

(2.20) (V/ + 7(v3VM)'',Vq) =0 yqeH\n), 

(2.21) {lpf,ly)+e{V^l^.yly) + {{ip[Vp + -fipV^i])^,Viy) ^0 e H\n), 

(2.22) (/»-£(v/,V^)- J((/(^))^^) =0 VyjeH\n). 

Setting ip = -ipf in ( |2.22| ), ^ fi^ in ( |2.2ip and = ^ in ( |2.20[ ) and adding the 
resulting equations we get 

Jt' ■ — - 7 

where 



(2.23) ^ J.(¥'') + £ II v/IL. + - ||Vp* + 7(^VM)i^. = 7^^(^) 



1 



n\t) := - (/(/) - ifiip))', ^f) + (V/ + 7(^Vm)^ V'V/ - (^Vm)O 
+ (^[Vp^ + 7(^Vm)*] - (v5[Vp + 7'/'V//])^ V//*) . 



Integrating (2.23) in t gives 

Mv\s)) + f^(e ||v/(0|r^. + ^ ||v/(t) + 7(^V/i)^(t)||^,) 



Jo 



Note that, for each fixed e > 0, f{^) e L2((0, T); i?i(0)), since e L''((0, T); L°°(r2)), 
and /( • ) is continuous. Sending i5 — > 0+ and using properties of the Steklov average 
Si (cf. 19, Ch. 2]) we get 



lim / n^{t)dt = 0, 
^0+ Ja 



|V/i(i)lli2 + - ||Vp(t) + 7(^(t)V/.(t))||2, ) dt = Je(^(0)). 



7 



Hence, we recover (2.121. Observe that from (2.12) we can conclude that E{t) is an 
absolutely continuous function of time. 
Using (2.12) and the estimate 

(2.24) (F((p),l) > ^llQjIlf. - 



(^^(^),l)>^ll^lli^ 



we discover that 



2 < C, for all time and independent of e, and inequality (2.131 



follows. Inequalities (2.14) and ( 2.15| ) follow straightforwardly from (2.10), (2.13), 
and the Sobolev embedding H^{n) ^ L^{n) for d = 2, 3. Inequality (2.161 is an 
immediate consequence of (2.9) and the fact that e || V/i||^2(i2) < C < oo from (2.12). 
To sh ow ( |2.17| ), by ^3^, the Sobolev embedding H\n) ^ L^{n) for d = 2, 3, and 
( plsl , we get for any ly e W^^^{n) 

{ipt.v) = -e(VAi, Vv) + {ifVL, Vv) 

< e||VM||^2 ||Vi/||i2 + 11(^1^6 ||U||^2 ||Vi.||i3 

< C[e\\V^l\\L2 + M\Hl \MlA I|Vi^||l3 
<C e||VAi|U2 ^ 
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The above inequality and (2.12) infer (2.17) 



Finally, (2.18) is an immediate consequence of (2.12) and the inequality {(p — 
1)^ > {W\ — l) ■ The proof is complete. □ 

The next lemma shows that weak solutions have some additional regularities. 
Lemma 2.3. Suppose that e H^{n) and Vl C M.'^ {d = 2, 3) is a Lips- 
chitz domain. Let (p, /i, 1^9) be a weak solution defined by (2.8) -(2. 10). Then ip G 



((0, T) ; (ri)) . Moreover, if fl is a convex polygonal or polyhedral domain, then 



Proof. We begin by rewriting (2.10) as 



(2.25) 



Hence, is a weak solution to a Poisson equation with homogeneous Neumann bound- 
ary conditions and the right-hand side "source" function g :— fi ~ f(ip). Since 
if e ((0,T) (rj)), then 5 G ((0,r) (O)). By elliptic regularity theory 
(cf. [TO, Ch. 7]) we conclude that Lp {{0,T);H^ (fl)) and = on dn in the 
distributional sense. 

Introduce the function space 

09 



t/ := |6I e H^{n); ^ = on dn} 



For any 61 e C/ n C°°{n), setting ^ = 
density argument we get 

e{Aip,Ae) = (V [^i- 



which implies that (p e H^{il) with {p, 1 
following biharmonic problem: 



A9 in (2.25), integrating by parts and using a 



-^f{p)],we) yeeu, 

= const is the unique weak solution to the 



dp 
dn 



AV = -Ag 
dAp 
dn 



= 



in il, 
on dfl. 



Since —Ag E ((0,T) ;H^^ (^))i it follows from a well-known regularity result of 
[5] that p € L"^ ((0, T) ; (fi)) when £7 is a convex polygonal domain. The proof is 
complete. □ 

We conclude this section by establishing the following uniqueness theorem for 



weak solutions of problem (2.1), (2.2), and ( |1.6[ )-( 1.8 1 defined in Definition 2.1 

Theorem 2.4. Suppose that p^ £ H^{fl) with Js (Po) < Co for some e- 
independent constant Cq and C M'* {d — 2,3) is a Lipschitz domain. We say that 
a weak solution (p, /i, p) belongs to the function .space T if it satisfies the additional 
regularity conditions Vp + -fpV^jL e L^{{Q,T);L'^{^1)), fi e {{0,T); (Q)) , 



and pt & L^((0,T); {H^{Q))*). Then weak solutions of (2.8)-( 2J^ in the function 
class J- are unique. 

Proof. Since the proof is long, we divide it into six steps. 

Step 1: Suppose {pi, fii, pi), i = 1,2, are two weak solutions, and define u,; := 
-Vpi~jPiVfii, i = I, 2. Letp = pi~P2, u = U1-U2, ^ = fii-fJ.2, and , p = px-p2- 



Subtracting the corresponding equations of (2.8)-(2.10) satisfied by (j)\,[i\,p\) and 
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{P2i fJ'2i we get the following "error" equations: 

(2.26) (u,Vg)=0 e H^n), 

(2.27) {ipt,iy)+e{V^i,Viy) -{ipiu + ipu2,Viy)^0 ^ly e H\n), 



(2.28) 



We will frequently use the fact that f{x, t) dx — 0, for almost all t e (0, T), which 
follows from (2.111. 

Setting 1/ = in ( 2.27[ ) and -0 = in (2.28), adding the resulting equations, and 
using the fact that (u2, V (<y3^)) = = (u, \/{Lpiip)) we get 

\jt " -(u- + ^(5((Pi,V52)V,m), 

where g{(pi,ip2) ^^p\ + 'Pi'P2 + V?i - 1. 

Using Schwarz inequality and the following Gagliardo-Nirenberg inequality (cf. 



in (2.291 we get 



^ + 2 ||/i||^2 < 2||u||£2||V^i||L2||(^||i=o + ^||g(v5i,(^2)||L- ||¥'|Il2 ||Ai|lL2 

e C 2 

< ^I|u||i2 + -||V(pi||i2||A(^||22||<^||J + -\\g{^u^2)\\L^ M\l-^ IImIIl2 



£ 

47 ' 



< -l|u||i2 + |^||A^||i2+C(e)||V^i||-' 



2 

L2 



+ ^ll5('^l,V'2)||ioo \\ip\\\2 



Hence, it follows from (2.12) that 
d 



(2.30) 



dt 



L2 



< 



U 



47" ■ 16 

+ C(e)(l + 115(^1,^2)11 



2 

L2 



ll^ll 



L2 • 



Step 2: Setting ^p = Aip in ( |2.28| gives 



Hence 
(2.31) 



e||A(^lli2 



(m, A(^) + ^ (g((^i, (^2)<^, A(p) 
1„ „o 1 



llA^lli2<-iiMlli2 + c(s) 11.9(^1,^2)111.. Ilv'lli. 



3: First, on noting that 

U = Ui - U2 = -Vp - 7((piV^i - <y92V^2) = -Vp - ^iflV fi ~ -JipV fJ,2, 
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and using (2.26) with q ~ p we obtain 
(2.32) 



= - - ((^u, v^2)- 



Second, setting v = ji in (|2.27|) yields 
(2.33) 



(V3t,/i) +e\\Vp.\\l2 = ((y5lU + (^U2,V/i). 



Third, applying the Steklov average operator S\_ to ( |2.28| ) (we use the same notion 
as in the proof of Lemma 2.2) yields 



(m^^) -e(V(p^v^) - - ((5(^1, ^2)v)',^) =0 

Where 17(1^1, (^2) ~ ^1 + + </'2 ^ 1- Setting = —ipf in the above equation gives 

- + + J {{g{vu^2)v)\vl) = 0. 

Taking the limit (5—5-0+ and using the properties of Steklov average operator (cf. 
P] ) we get 



(2.34) 



e d 1 
z at £ 



Finally, adding ( |2.32D , ( |2.33[ ), and ( |2.34D , using the fact that (u2,V(<p/z)) = 
(u, V((^^2)) — and Young's and Gagliardo-Nirenberg inequalities (cf. [TJIH]) we 

get 



= ((^U2, V^) - (V3U, V^2) = -(U2 • VyJ,^) + (u • V(y5, ^2) 
< ||u2||L2||V'^||L3||/i||L« + ||u||L2||V</'||L3||Ai2||L6 



l|V<^||i2 + 7 ((^t,g((^i,v32)</?) 



< 



C{\\n2\\LAMH^ + ||u|U.||/i2||ffi)(||A^|||,||V¥>||2^'* + llVy^ll 



<-llul|i. + ^||A*r 



47 



16 



IIA^II 



L2 



C(£) ||u2|12-+||a.2||^-/ liv^lli. 



Hence 
(2.35) 



47 



u 



ii. + ^iiv/.|ii. + ||iiv^iii. 



((^f,5r(<^l,(^2)<^) 



<^||A^||i. + |||M||i.+C(£)(||u2||-^+||M2|| 



livy^lli.. 



S'tep ^: To control the last term on the left-hand side of (2.35), we rewrite 

5(<^l, (^2) + ifil - I = {(pi ~ ip2f + ipiV2 - 1 ((^5^ - 1) + i<PW2- 
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Thus 



1 3 
(2.36) {ipt,g{ipi,ip2)(p) = 2 {{v'^)t, V'^ - l) + 2 (('/'^)t' V'l'^a) 



d 

dt 
d 

dt 



\\Mh + ^ (<p^</'l<P2) 



(('^l</'2)t,'/?^) 



To bound the last term on the right-hand side of (2.361, we mtroduce the inverse 



Laplace operator A'^ : {H^{n)y H^{n). For any w e iH\n))* with (it;, 1) = 0, 
let A~^w G H^{n) be the unique solution of the following problem: 



(2.37) 
(2.38) 



(A~^w,l) = 0. 



It is straightforward to show that, for all w e {H^{n))* with {w, 1) = 0, 

(2.39) \\w\\^H,y^\\V{A-'w)\\^,. 

Then for j,k = 1,2 and j ^ k using Sobolev inequality (of. T]) we have 

(2.40) {^jt, fkv'') = - (V(A- Vjt), '^'V^fc + 2^k^V^) 

< ||V(A-Vjt)IUHI<y5^V¥'fc + 2^fe(^V^||L2 

< \\^jt\\(^my (l|V(y5fe||L6||^||^6 + 2||^fc||L6||(^||i6||V(^||L6) 

< C\\ipjt\\^H^y (^WAipkWL-Al^fWh + Wy^kWn^ l|V<y9|U2||A(p||i2^ 



< 



-\\Aip\\l2 +C{e)(\\ip.jt\\^H^y ||A95fc||i2 + WipjiWln^y llv'fellHi) ||V(/3|| 



2 

L2- 



step 5: Adding ( |2.30[ ), ( |2.31[ ) and e times of ( |2.35[ ), and utilizing ( |2.36[ ) and ( |2.40 | 

we get 



(2.41) 



1 13 

2 II'^IIl^ + yl|V(^||i2 + -||<^||i4 + - {ip^,<fiiip2) 

+ 2(1 - llA^lli^ + 1^ IIHli. + y llVA^IIi^ + f^ll A^||i2 



<C(e)(l + ||5((^i,^2)|| 



C(e) 



|u2!i£r' 



1^^211 



j.k = l 



<a(i) (11^11^2 + ||V(/.|| 

L2 ) I 
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where 



|u2||l^'^+||M2||^r 



a{t):=C{s) l + ||.g((pi,^2)IL. 



j,k=l 



Integrating (2.41) in t over the interval (0,t) we get 



(2.42) 



II^WII 



L2 



- e'll v^(t)||i. + Mt)\\l, + 6 {^\t), Mt)Mt)) 



< 



dt. 



Step 6: Define 



max 



[te[0,T]; {ip^is),ip,{.s)ip2is)) >Oys e [0,1]}. 



We now show that t > 0. Since (pi{0) — 'P2i0) — ^Pq, by continuity there exists ti > 
such that for j = 1, 2 



1 



a2' 

Consequently, 



(2.43) 



{ip''it),ipiit)p2is)) > IIMI 



,e||2 



L2 



12 



12 



,ri] yte[0,ti], VO < 77 e L°°(17). 



MI2 yte[0,ti] 



Substituting ( |2.43| ) into ( |2.42| ) yields 
1 



(2.44) 2"^^*)"^^ 



< 



2 , iivv../„Ml2 



«(s)(ii^G^)irL^ + iiv^(s)iii 

By Gronwall's inequality we get 



L2 

ds 



yte [0,ti]. 



(2.45) 
1 



Mt)\\l.+e^W^ml,.< 



11^(0)11^. +£2||V^(0)|| 



L2 



cxp ■ 



a{s) ds 



Vie[0,ti]. 



Here we have used the fact that a{s) ds < 00. Thus, tp{t) = for f G [0,ii]. 
Therefore, t > ti > 0. In fact, the above proof also shows that (f{t) = for t € [0, t]. 

Suppose that r < T, by the definition of t we have <^(t) = 0, that is, — 
(P2{t). Repeating the above Gronwall's inequality argument with r in place of t = 0, 
we conclude that there exists t2 > t such that (p{t) — ioi t E [0,^2]- Hence, 
(^(p'^{s),(pi{s)(p2{s)) = Vs G [0,^2]- By the definition of r we must have r > ^2. 
So we get a contradiction. Therefore, r = T and ip{t) — 0, i.e., ipi{t) — ip2it), for 
t e [0,T]. The proof is complete. □ 
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3. Fully discrete finite element method. 



3.1. Formulation of the finite element method. For simplicity we assume 
that Q C R'^ {d = 2, 3) is a polygonal or polyhedral domain. Let ~ {t„i}m=o ^ 
uniform partition of [0,T] of mesh size r := and dtv^ := (tj™ — u™^^)/t. (This 
is for simplicity; we could also use a quasi-uniform partition of the time interface.) 
Let be a quasi-uniform "triangulation" of the domain Vl of mesh size h e (0,1) 
and n = UifGTh ^ ^ tetrahedrons in the case d = 3). For a nonnegative 

integer r, let Pr{K) denote the space of polynomials of degree less than or equal to r 
on K, and define 

SI = {vh e C%n) I Vh\K e Pr{K) VK e %] ■ 

For fixed positive integers r and we introduce the finite element spaces Vh = 3"^^ 

and Wh = Sh- Define Vh ■— Vh n Lo(^) ^^'^ similarly for Wh- 

We now are ready to introduce our fully discrete finite element method for prob- 



lem pll), (2.2), (1.6)-(1.8) based on the variational formulation (2.8h(2.10|. Find 



{ipT,IJ'h,<T}^=i CWhX Vh X Vh such that 



(3.1) 
(3.2) 

(3.3) 
(3.4) 



idt^^,:.h)+e{V^i'^,V,.h) 

+ {^h~'[Wh + T^r^v^r], Vi.,) = 

ifiT,^Ph)-e{V^^,V4'h) - - (/r^V-zO = 



yqh e Wh, 

yi^h e Vh, 
y^Ph e Vh, 



VOh, 



where ipoh G Vh, to be specified in the next section, is an approximation of (pg, and 



(3.5) 



ih ■= Wh ) 



Vh ■ 



We will prove that {Wh, V^,, T4) is a stable triple for our mixed finite element approxi- 
mation. The techniques we use are based on energy estimates and convexity analysis, 
rather than an inf-sup-type condition, which would be used for the analysis of linear 
biharmonic-type equations [9l [TOl [TH |26] . 



3.2. Well-posedness of the finite element method. The goal of this sub- 



section is to show that the fully discrete finite element scheme (3.1|-(3.5| is uniquely 



solvable and energy stable for all h,T,e > 0. To prove unconditional unique solvability, 
we shall show that at each time step the scheme can be reformulated as a minimiza- 
tion problem for a strictly convex and coercive functional. We begin by defining an 
inner product on the subspace Vh- 

Lemma 3.1. Define the bilinear form a : x V/j M via 



(3.6) 



where Ai{e,^p) — e + ^Lp^ and p{^) G Wh solves 
(3.7) (Vp(m), Vg)^. = -7 V^., Vg)^, 



yqeWh - 



Then a{ ■ , ■ ) is an inner product on Vh- 

The proof is omitted for brevity. Note that in the next few calculations, the 
pressure, p, will be regarded as an auxiliary variable that can be calculated when the 
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chemical potential, ji, is known. Owing to the last result, we can define an invertible 
linear operator jC : Vh ^ Vh via the following problem: given ( gV^j find /U € V/j such 
that 

(3.8) a{^i,u) = -{C,u)^, \/i^€Vh. 

This clearly has a unique solution because a(- , • ) is an inner product on Vh- We 
write = — C, or, equivalcntly, ^ = —C~^{(). 

We now wish to define a negative norm, i.e., a discrete analogue to the norm. 
Again we omit the details for brevity. 

Lemma 3.2. Let (, &Vh and suppose fi^, /xj G Vh are the unique weak solutions 
to C (/K^) = — C and C (/i^) = — ^. Define 

(3.9) (COr-i :=a(MC.M£) = -(C,M«)i,2 = -(Mc>0i,2. 
( • , • defines an inner product on Vh, and the induced norm is 

(3-10) llCllr-i = ^/{COc-^- 

Using this last norm wo can define a variational problem closely related to our 
fully discrete scheme. 

Lemma 3.3. Set Ki := l)^,, and define tp™"^ := (^^"^ - Ki & Vh. For 

all (fi € Vh, define the nonlinear functional 

(3.11) G{^) := I 11^- + 1 y + K^Wl, + I llV^lli. - ^ (^r',^)^. • 

G is strictly convex and coercive on the linear suhspace Vh. Consequently, G has a 
unique minimizer, call it e Vh- Moreover, e Vh is the unique minimizer of G 
if and only if it is the unique solution to 

(3.12) \ («'+i^i)\v)^^ +e(V(^r,v^)^. - (Mr,V')i. = \ (v'r'.^)^. 

for all tjj GVh, where G Vh is the unique solution to 

(3.13) a «, 1^) = - - Vi/ e Vh. 

Proof. In detail, the first variation, i.e. the gradient, of the first term of G is 

= -(/^>V')i2 , 

where, owing to the definition of the inner product ( • , • )£-i, /x € is the unique 

solution to 

(3.14) a{f,,iy) = -{^-^T-\'^)L^ Vz. e I4. 



ds 
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The second variation is 



which estabhshes the strict convexity of the term. The strict convexity of G foUows 
because each of the other terms is at least convex. The coercivity of G foUows from 
an estimate of the form 



(3.15) 



where < Ci(e),C2(e) < oo are constants. By the standard theory of convex opti- 
mization, G has a unique (bounded) minimizer in Vu, call it 95™. Moreover, Lp^ is the 
unique minimizer of G if and only if it is the unique solution to A^pm (jp) = 0, for all 

ip G Vh, where 



A^it^) ^Giifi + sij) = -((^ + Xi)3,7A) +e(V^,V^)^3 

CIS s—0 £ 



(3.16) 



1 



The rest of the details follow from the definition of the inner product {■,-)c-^- D 

Finally, we are in the position to prove the unconditional unique solvability of our 
scheme. 

Theorem 3.4. The scheme (3.1)-(3.3) is uniquely solvable for any mesh param- 
eters T and h and for any phase parameter e. 



Proof. First it is clear that a necessary condition for solvability of (3.2 1 is that 
(3.17) (^^,1)^. = «-\l)^, =:ifi. 



as can be found by taking = 1 in ( |3.2[ ). Now, let (^i™, </3™) G x V/j be a solution 
of ([3l2l)-([3l3l). Define ip"^ ipT + Ki/\n\. Set 



(3.18) 



7n\3 



1 



and define /x™ := /i™ + A'2/|fi|. Then it is straightforward to show that (/i™,(p™) € 
Vh X V;i is a solution to (3.2|-(3.3). In fact, there is a one-to-one correspondence of 
the respective solution sets. Namely, if (Mr'V™) G ^/i ^ is a solut io n to (3.1 )- 
(1331) th en {fi ]^ - K2/M,'pT ^ ^i/M) e 14 x 14 is a solution to |3.12|-|3.13| . But 
(3.12)-(3.13) admits only a unique solution, which proves that (3.1 )-(3.3) is uniquely 
solvable. □ 

We now establish a discrete energy law for the numerical scheme that mimics the 
continuous version (2.12). 



Lemma 3.5. Let (p™, /x™, </5™) denote the unique solution of the scheme (3.2)- 



(3.5) and define uJJ 



(3.19) 



El 



Vp™ -7(^™"'VAi™; then there holds 
1 „ 



m— 1 



7 



2 

h llL2 



T 

4e 



2s' WdtV^Uh 



+ 2M^dt^nl+2\\dt^nl^ 
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for allO<i< M. Here E'^^ Je{^h) and ■ ) is defined m (Ol 



Proof. The desired estimate (3.19) follows from setting — in (3.1 ), i^h = 
in (3.2 1, Tph = —dt^p'f^ in (3.3), adding the resulting equations, using the identities 



[dt\\wvul^ + T\\dtw^nl^], 



m\2 



Il2 



m^2|l2 
L2 



+ 2y^dt^nh+2\\dt^nh], 



and applying the operator t J2m=i ^'^ combined equation. □ 

The discrete energy law immediately implies the following uniform (in e, h, r) a 
priori estimates for (p™, /j,™, (/?™). 



Lemma 3.6. Let (p^ , fi^ , tp"^) be the unique solution of (3.1)-(3.5) and define 
u;^' := -Vp^" - "fifl'^'^V . Suppose that E° < oo. Then, for all m > 1, 



(3.20) 



(fldx, 



and, in addition, there hold the following estimates: 



(3.21) 
(3.22) 

(3.23) 
(3.24) 
(3.25) 
(3.26) 
(3.27) 



max 

0<m<A/ 



e\\Vipnl^ + -{F{^^),l) 



e 



< C. 



max < 

0<m<M £ 



M 

-E 

m— 1 



1 



M 

E 



M 

\\\'PhWh-Vh )\\L2 + Wh) -Wh ) 



IL2 



m— 1 



< c, 



< c. 



< Ce, 



c 



-E llVpnili 



A/ 
m— 1 



/or some e, h, and t -independent constant C — C{Ef^) > 0. 

Proof. (3.20) follows immediately from setting i^^ = 1 in (3.2), and ( |3.21[ )"( [3.25| ) 
are the immediate corollaries of the discrete energy law (3.19). (3.26) follows from 
the identity VpJ^* = -u™ - 7(p™- V^i^J", th e Sob olev embedding H^) ^ L^ifl) for 
d = 2, 3, Young's inequality and estimates (3.22) and (3.23). 

Let Qfi denote the standard projection operator into Vh (cf. [71 [3]), and for any 
1/ (z VF^'^(ri), set I'll ~ QhV in (3.2). To prove (3.27), we use the Schwarz inequality 
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and the Sobolev embedding H^{n) ^ L^{Q) (for d = 2, 3) to we get 



<£\\^^^h\\L-\\^Qh'y\\L■^ + \\^h 



m— 1 1 



< c 



h llL2 " 

eiivA^riL. 



1 



|VQhi^lli3 



L3 ' 



where we have used the W^'^ stabihty of the projection Qh (cf. [71 [9]) to get 
the last inequality. (3.27) now follows immediately from the above inequality and 



estimates (3.21) and (3.23). The proof is complete. □ 

Remark 3.1. Property ( 3.20[ ) says that the proposed numerical method enjoys 
the same mass conservation law as the phase field model (cf. (2.11 )j. This property 
will be validated numerically in Sec. 

4. Convergence analysis. The goal of this section is to prove that the fully 
discrete finite element solution has a unique accumulation point (in some function 
space) and this accumulation point is necessarily a weak solution to problem (2.8)- 
(2.10). A byproduct of this convergence result is to provide a constructive proof of 



the existence of weak solutions to problem (2.8)-(2.10) 



First, we derive some additional estimates for the finite element solution. To the 
end, we introduce the discrete Laplacian A;, : Vh which is defined as follows: 

for any Vh €Vh, ^h^h G Vh denotes the unique solution to the problem 

(4.1) {AhVh,Wh) ^ -{'^Vh,\7wh) VwheVh. 

In particular, setting Wh = A^tj/j in (|4.1[), we obtain 



\A.hVh\ 



L2 



(Vvh.VAhVh) 



Lemma 4.1. Let (p^ , fi^^ , (p^) and u™ be same as in Lemma 3.6 Then, under 
the assumption E^^ < oo, there hold the following additional estimates: for d = 2,3, 



M 



(4.2) 
(4.3) 
(4.4) 
(4.5) 
(4.6) 



Tj2\\^hVm^<Co{T+l)Cie), 



M 



rY^yTh-' <Co{T+l)C{s), 

m— 1 
M 

rJ2\\W^'nl<Co{T+l)C{e), 

TEl|vK'llf^<Co(r + i)c(£), 

— 1 

^El|dt^ril(1lf <Co(r + i)C(£), 



for some e,h, and t -independent constant Cq = Co(-E)J) > and some h and t- 
independent constant C{e) > that grows like e^^ , for some r G Z+ , as e — > 0. 
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Proof. Setting — ^h'-Ph^ in (3.3), using the definition of /S.h^™ , and the 
Schwarz inequahty we get 



Hence, 
(4.7) 



1 fe 



^ ii/riiL2 



<l||V^"'"^ 



?n 1 1 
h II L2 



ei!A,^5riii2<iivA^riii^ 



h llL2 



iii/nii. 



To bound ||//^" ||/^2, we write 

-f^Ti . /,„m\3 ,^m— 1 , ( ( , jm\2 -i \ i , ,„m 



Then by (3.211 wc have 



W £m\\2 ^oIL ^^'^ II 2 II / ,„m\2 i || 2 , || „m , „m— 1 ii 2 



L2 



= 8|i^riii== (F(^r),i) + ii¥'r"<"'iii= 
<ce\wni^ + \wt-vi 



m-l||2 

Ll2 



We now appeal to the following discrete Gagliardo-Nirenberg inequality (cf. ^ and 

mmy- 

d 3(4-ei) 

(4.8) M:U^<C\\/^HvXt'' \WnW^' + CWh\\l^ (d = 2,3) 

and get 

(4.9) 

d 3(4-d) 



ry^W, ^711 1 1 2 I II , , „(. 



m-l||2 

Il2 



< e 



< 



-\\A,^T\\h + Ce^^ + Ce WnWle + hi 



11^ 



II A.^nii^ + c £^ + £ ii^riii. + ikr - 



lL2 



lL2 



where we used the Sobolev embedding H^{n) ^ L^(ri) for d = 2, 3 in the last step. 
The n (|4.2[ ) follows from applying the operator r ^^^-^ to ( |4.7[) an d using (4.9), (3.22), 
and (3.231. (4.3) is an immediate consequence of (4.8 1 and (4.2). 

To prove (4.4), we recall another discrete Gagliardo-Nirenberg inequality (cf. [TH] 
and Din]): 



(4.10) WVi^hh. <C\\Viyh\\J \\Ahiyh\\l2+C\\Vi^hh2 Vj/,, e y,„ d = 2, 3. 
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It follows from the above inequality and estimates (3.211 and (4.2) that 

M 



M 



2(i-d) 



M 



M 



(4.11) 



<^rJ2\\^HVH\\h + ^T 

£ d £d 



which proves (4.4 1 



Inequality (4.5) follows from the estimate 



4(6-d) 4(6-d) 

,m 1 1 12-d [- I 



L2 



4(6-d) 4(6-d) 
m-l|| 12-d ||Y7,,m|| i2-<i 
'Ph IIloc II VAt/i 11^2 



< Cllu 



h llL2 



7n — l I 



4(6-d) 



and estimates (3.23) and (4.3) 



Now, let Qh denote the standard projection operator into Vt (cf. [HE])- For 
any v G H^{il), setting i^/j = Q^i/ in (3.2), we get 



< 



£iiv/iriu= 



I , ^m—l II II II 



< c 



e\\v^^^nL^ + y^r'\\L-\\<\\L■^ 



llVi^lU^, 



where we have used the stabihty of the projection Qh (cf. [3E|) to get the last 
inequality. (4.6) now follows immediately from the above inequality and estimates 



(4.3) and (3.23). The proof is complete. □ 



Next, let (pfi.T{x,t) denote the piecewise linear interpolant (in t) of the fully dis- 
crete solution that is. 



(4.12) 



y^Kri- ,t) : = 



t - t 



m— 1 



for 1 < m < M. Let p^ .^{x,t), Uh,T{x,t), fif^ .^{x,t), (pf^ .^{x,t), and ip^ .^{x,t), de- 
note the piecewise constant extensions of {p™}, {u™}, {/i}"}, and {i^a™}, respectively, 
defined as follows 



(4.13) 




71™ 

Ph 




- 1 : im] 7 


1 < m < M, 


(4.14) 


Uft,r(-,i) 


— 11™ 




- 1 : im] 7 


1 < m < Af, 


(4.15) 




— //™ 




- 1 1 im] 7 


1 < m < M, 


(4.16) 




:= 'Ph 




- 1 1 im] 7 


1 < m < M, 


(4.17) 




,^m — 1 




-1 : ^m] 7 


1 < m < Af. 



We remark that Lph.rix, t) is a continuous piecewise polynomial function in space and 
time, py^^{x,t), u/i.T-(a;, t), JIf^_^{x,t), and Jpf^ .^{x,t) are right continuous at the nodes 
{tfn}, and Tpi^ ^{x, t) is left continuous at the nodes {^m}- 
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The main result of this section is the following convergence theorem. 
Theorem 4.2. Let ft C M'' (d — 2,3) be a polygonal or polyhedral domain. For 
each fixed e > 0, suppose that J^dfoh) l£ Co < oo, where Cq is independent of h, and 

lim llm - flh^ = 0. 

Then the sequence {(p^ u^^r, M/i r i ^/i r)} ^^■^ '^'^ accumulation point (p"^ , , /i*^ , ) 
with = — Vp^ — ^if^V , and {p^ , ^'^ , ip"^) is a weak solution to problem (2.8)- 
( [2T0l ). 

Proof. We divide the proof into two steps. 

Step 1: Extracting convergent subsequences. The estimates of Lemmas |3.6| and 



4.1 immediately give the following (uniform in h and r) estimates: 





^(L2) + 


l^lr-l| 














<c, 






+ W'^Ph.r 




<c, 


















La (L") 




L2((M/i 


an + l 


i^h,T)t\\L 


'((ffi)*) 





(4.18) 
(4.19) 
(4.20) 

(4.21) 
(4.22) 

(4.23) II (9.;, 

where /3 :— '^^^^'^'^ > 4 and a := ^j^^^-* > |. Note, we have suppressed the de- 
pendences of the constants on T and e above. {Vp^ ^-l is uniformly (with respect 
to h and r) integrable in L"' ((0, T); L^(J7)) and {(<y5/i,r)t} is uniformly integrable in 
L'^ {iO,T);{H^r). _ _ _ 

Then there exists a convergent subsequence of {(p^^ ^, U/j^t-, (p/j^,-)} (still de- 
noted by the same symbols) and a quadruple (p^, u*^, /i^, (/s^) such that 

p' e {io,T);H\n)nLlin)) , e ((o,r);L2(r!)) , 

G ((0,T);i/i(r!)) ni4'^((0,r);L-(r!)), e ((0,T);Tyi'4(^)) ^ 

e ((0, T); (i/i(f^))*) n ((0, T); (l¥i^3(j^))*) ^ ^ ^2 ^^q^ T); H^Q)) , 



' weakly in L" ((0, T);H\n) D Ll{n)) , 

weakly in {{0,T);L^{n)) , 

weakly in {{0,T); H^{n)) , 

(fi' weakly* in L°° ((0, T);H\n)) D ((0, T); i°°(n)) , 
strongly in {{0,T); L^{n)) , 

weakly in ((0,T); (M^^'^^^))*^ ^ p^i^a ^^q^ j.^. (iyi(i7))* 
weakly in {{0,T);W^'^{n)) . 



and 






(4.24) 


Ph^r 




(4.25) 


Uh,T 




(4.26) 


M/i,r 




(4.27) 







We have used Aubin-Lions lemma (cf. ^5^) to conclude (4.27). 
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From (3.24) we also have 



M 



\'fh -fh \\m 



ni—l 
M 



T 



dt 



\ ^ \\,„ra ,„m-l||^ 

o 2^ ll<^/i -Vh \\m 



1||2 T\,p 



Hence, {Lph,T}, {^h r}; ^^^d {(^^ ^} converge to the same Umit as h,T ^ 0. 

Step 2: Passing to the limit. We now want to pass to the hmit in (3.1 1-(3.4| and 
to show that {p'^ , , (p'^) is a weak solution to problem (2.8)-(2.10) with the initial 
data </3^(0) — ifQ- To this end, we rewrite (3.1)-(3.4) as 

(4.28) {uh,r,'^qh) =0 yqh^Wh, 

(4.29) ii^h,T)t, Vh) + e(V/i/i,^, Vi/,,) - (p;,,^u,,^^, Vvn) = Vi/^ e V,,, 



(4.30) 



where denotes the right continuous constant extension of {/™}. 

For any 77 G C°([0,r]), multiplying (4.28)-(4.30) by ij, respectively, and integrat- 
ing the resulting equations in t from to T we get 



JQ 

-{^f, .^Uh,T,'^Vh)^vit)dt = e Vh, 



(4.31) 
(4.32) 

(4.33) 



For any {q,iy,tP) G [H\n) D C\n)]^, let {qh,i^h,^h) € Wh x Vh x Vh be the 
standard finite element (nodal) interpolations of {q,v,ilj) in ( |4.31 1-(4.33 1. Since 

qh — > q. Vh — > V, iph — > V strongly \n H (il), 



sending /i, r ^ in (4.31 )-( 4.33 1 and using (4.25|-(4.26l we get ip^{Q) = ipl and 



/ (u% Vg)?7(t)rfi = VqeH\n), 
Jo 



(4.34) 

(4.35) |((^^,^^)+£(V^^V^/) + (v3''u^Vl')}?7(^)rf^ = Viy e H\n), 

(4.36) ^^{(/i^V')-£(V^^VV')- J(/(^^),V)}?](t)rfi = y^j&H'in). 
Moreover, from the identity u™ = —Vp^h ^ IvT^^'^ I^T '^^ have 



(uh,r,'^q)v{t)d.t 



f 

Jo 

/ i:^h,r»-h,T,^v)ri{t)dt^ - I (Tpf^,^[Vpf^.^ 
Jo Jo 



{^Ph,r + l'Ph,T^l^h.r-,^q)ri{t) dt, 
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Sending h,T ^ and using (4.24)-(4.26) yields 



(4.37) / (u^Vg)r;(^)d^-- / {Wp' + j^'Vfi' ,Vq)f]{t) dt, 
Jq Jq 

(4.38) [ {if''\f.yv)ri{t)dt^-( {ip^l^Ip^ +-iip''\7ii%\7v)ri{t)dt. 
Jo Jo 



Combining ( |4!34l )-( |4!36l ) and ( [437| -( |4!38l ) we obtain ( [zl|)-([2lol), si nce C°[Q,T] 
is dense in L^(0,T). Hence, (p^,/i^,(^^) is a weak solution to (2.8)-(2.10). The proof 
is complete. □ 

Corollary 4.3. The whole sequence {{pf,_^,Uh^Ti'Ph r^Vh t)} converges if weak 
solutions to problem (2.8)-(2.10) are unique. 



Proof. We have shown in the above proof that {(p^ ^, fihr,ipf^^}ha.s a convergent 

-( |2.10| ). Moreover, the 



subsequence and its limit {p^ , (^^) is a weak solution of (2 



proof also implies that the limit of every convergent subsequence of {(p/j ^, /I^ .^i t)} 
is necessarily a weak solution of ( 2.8 H( |2J^ ). Hence, by the uniqueness assumption 
of weak solutions, the whole sequence {(p^^, ^, (/s^^)} must converge to the unique 
weak solution. □ 

Remark 4.1. In Theorem 



4.2 and Corollary 4.3, ft is assumed to be a polygonal 



or polyhedral domain. This assumption is imposed only to avoid the technicalities 
for defining our finite element method (3.1 )-([3^. It is not used or needed in the 



proofs of the theorem and the corollary. By using the standard numerical integration 
technique or the approximated boundary technique (i.e., to approximate a bounded 
Lipschitz domain by a sequence of polygonal or polyhedral domains) (cf. it can 

be proved that the modified finite element methods would also possess all the properties 
proved in Lemmas \3.1\ \3.S\ \3.S\ \3.^ \3.(^ and Theorem \ 3.4\ as well as Lemma \4.1\ 
As a result, the conclusions of Theorem \4.S\ and Corollary \4.3\ still hold when fl is a 
bounded Lipschitz domain. 

From Theorem |4 . 2 1 and Theorem |2.4| we immediately have 



Theorem 4.4. There exists a weak solution to problem (2.8)-(2.10) and weak 
solutions are unique in the function class T . 

We conclude this section with a remark on the error estimates for the solution of 



the fully discrete scheme (3.1|-(3.4|. Using the standard (perturbation) technique as 



presented in |13j . it is not hard to prove that the scheme converges optimally in the 
energy norm. However, the error constant would contain a factor of exp (e"^). Such 
an error bound is clearly not very useful for small e. A better error bound would 
only depend on in some low polynomial order (cf. |14[ ll5jV Deriving such a 
polynomial order error bound is an on-going project and the result will be reported 
in a forthcoming paper. 

5. Numerical experiments. In this section we provide some numerical experi- 
ments to gauge the accuracy and reliability of the fully discrete finite element method 
developed in the previous sections. For the experiments we take = — Sj^ for 
simplicity. We use a square domain ft — (0, 1)^ and take 7^, to be a regular trian- 
gulation of $7 consisting of right isosceles triangles, as depicted in Fig. |A.1| We use 
a nonlinear multigrid method, which is detailed in Appendix lAl to solve the scheme 



(3.1)-(3.4| at each time step. We perform a battery of three tests on the scheme. 



First, we measure numerical convergence of the scheme in the presence of added, ar- 
tificial source terms. Second, we measure the numerical convergence of the scheme 
without source terms using a Cauchy-convergence method. Third, we conduct a test 
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h 






rate 


IIcmIIl^ 




rate 






rate 




8.683 X 10- 


3 




1.088 X 10- 


-2 




1.270 X IQ- 


-2 




V2/32 


1.850 X 10- 


3 


2.23 


2.701 X 10- 


-3 


2.01 


2.479 X IQ- 


-3 


2.35 


V2/64 


4.568 X 10- 


4 


2.01 


6.759 X 10- 


-4 


2.00 


5.759 X IQ- 


-4 


2.11 


^/l28 


1.141 X 10- 


4 


2.00 


1.691 X 10- 


-4 


2.00 


1.413 X IQ- 


-4 


2.03 


v^/256 


2.852 X 10- 


5 


2.00 


4.227 X 10- 


-5 


2.00 


3.515 X IQ- 


-5 


2.00 


Table 5.1 



convergence test. The final time is T = 1.0, and the refinement path is taken to be r = 
25.6?t^. The other parameters are e = 7 = 1.0; O = (0, 1)^. The global error at T is expected to be 
0{t) + O (h^) = O (h^), and this is confirmed. 



h 






rate 






rate 






rate 


V2/16 


2.886 X 10- 


1 




2.907 X 10- 


1 




2.943 X 10- 


1 




V2/32 


1.455 X 10- 


1 


0.99 


1.462 X 10- 


1 


0.99 


1.466 X 10- 


1 


1.01 


V2/64 


7.290 X 10- 


2 


1.00 


7.320 X 10- 


2 


1.00 


7.313 X 10- 


2 


1.00 


v^/l28 


3.647 X 10- 


2 


1.00 


3.660 X 10- 


2 


1.00 


3.653 X 10- 


2 


1.00 


^/256 


1.824 X 10- 


2 


1.00 


1.839 X 10- 


2 


1.00 


1.826 X 10- 


2 


1.00 


Table 5.2 



convergence test. The final time is T = 1.0, and the refinement path is taken to be t = 1.6h. 
The other parameters are e = 7 = 1.0; C = (0,1)^. The global error at T is expected to be 
0{t) + 0{h) = 0{h), and this is confirmed. 



of spinodal decomposition using varying values of the excess surface tension 7, and 
demonstrate the discrete energy dissipation and mass conservation properties of the 
scheme. 

For the convergence of the problem with source terms, we solve a problem of the 
foUowing form: find (p™, fT) ^ x Vh x Vh, such that 

(5.1) (Vpr + 7^r"'V/^",Vg^) = (si(x,i„),(j;,) V % G Vu, 

(5.2) {dt^l\i,n) + e{V^i'j:,Vvu) 

+ (^™-i[VK" +7¥'r~'V/ir], V;.^) - (,s2(x,i„),^.;,) V e Vh, 

(5.3) (Mr,V'/.)-£(V(pr,vv',.)- J(/r,V'h) = (s3(x,i„),v/^) vv^e^/., 

(5.4) = 

for m = l,...,Af, where the source terms are chosen so that the solution of the 
corresponding continuous problem is precisely 

(5.5) p{x, y, t) = fj,{x, y, t) = ip{x, y, t) = cos(7ri) • g{x) ■ g{y), 

with g{£^) = 16^^(^ — 1)^. The initial data are precisely given by (poh = {'■pi ' >0)): 
where Ih : {^) — > Vh is the standard nodal interpolation operator. All integrations 
are done exactly using the appropriate Gauss-quadrature rules. This is of course made 
possible since we are using polynomials in space. The exact values of all of the other 
parameters used in the test are given in the captions of Tabs. 5.1 and 5.2 Th e res ults 
of an error analysis using a quadratic refinement path are found in Tab. 5.1 and 
confirm the expected optimal second-order convergence rate in this case. The results 
of an error analysis using a linear refinement path are found in Tab. 5.2 and 
confirm the expected optimal first-order convergence rate for this case. Notice that 
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rate 


Pp|li2 




rate 


v^/l6 


^/32 


5.514 X 10- 


-2 




2.890 X IQ- 


1 




3.099 X 10- 


-2 




V2/32 


V2/64 


2.165 X 10- 


-2 


1.35 


1.229 X IQ- 


1 


1.23 


1.148 X 10" 


-2 


1.43 


V2/64 


^/l28 


6.284 X 10- 


-3 


1.78 


3.588 X IQ- 


2 


1.78 


3.250 X 10- 


-3 


1.82 


^/l28 


^/256 


1.636 X 10- 


-3 


1.94 


9.327 X 10- 


3 


1.94 


8.420 X 10- 


-4 


1.95 


v^/256 


v^/512 


4.132 X 10- 


-4 


1.99 


2.355 X 10- 


3 


1.99 


2.128 X 10- 


-4 


1.98 



Table 5.3 



Cauchy convergence test. The final time is T = 4.0 X 10 ^ , and the refinement path is taken 
to be T = 1.024/i2. The other parameters are e = 6.25 X lO'^; 7 = 1.25 X 10'^; U = (0, 1)^. The 
Cauchy difference is defined via 5^ := ipf^ — (ph^, where the approximations are evaluated at time 
t = T, and analogously for (5^, and 5p. The norm, of the Cauchy difference at T is expected to be 
0{r) + 0(/i,2) =0{h'^). 



he hf \\S^\\jji rate P^ll^i rate ||<5p||/fi rate 





V2 


/32 


8.569 X 10- 


1 




1.301 X 10- 







8.371 X 10- 


2 






V2 


/64 


4.160 X 10- 


1 


1.04 


6.295 X 10- 


1 


1.04 


3.715 X 10- 


1 


1.17 


v^/64 


\/2/ 


128 


2.061 X 10- 


1 


1.01 


3.111 X 10- 


1 


1.02 


1.779 X 10- 


2 


1.06 


v^/l28 


V2y 


256 


1.029 X 10- 


1 


1.00 


1.554 X 10" 


1 


1.00 


8.834 X 10- 


3 


1.01 


v^/256 


72/ 


512 


5.146 X 10- 


2 


1.00 


7.777 X 10- 


2 


1.00 


4.422 X 10- 


3 


1.00 



Table 5.4 

Cauchy convergence test. The final time is T = 4.0 X 10~^, and the refinement path is taken 
to be T = 2.0 X 10-^h. The other parameters are e = 6.25 X 10~^ ; 7 = 1.25 X 10'^; n = (0, 1)^. 
The norm of the Cauchy difference at T is expected to be 0{t) + O (h) = O (h). 



the approximations p™, (^™, and /i™ all appear to converge at the same optimal rates, 
in both cases. 

We now give the resuhs of a test without any artificial sources. In other words, 



we solve the scheme (5.1)~(5.3) with Si = 0, i = 1,2,3. The initial data are taken to 
be 

f [l.0-cos(4.0^x)] ■ [l.0-cos(2.0^y)] \ 
(5.6) ipoh = - 1-0 , 



and the parameters are given in the captions of Tabs. |5.3| and |5.4| Note that in this 
case we are not in possession of the exact solutions. To circumvent this, we measure 
the difference of the computed solutions at successive resolutions. Specifically, we 
compute the rate at which the Cauchy difference := ipf^^^ — ipff" converges to zero, 
where hf — 2hc, Tf = 2Ptc {p — 1 for a linear refinement path and p = 2 for a 
quadratic refinement path), and TfMf ~ TcMc = T. A quadratic refinement path, 
i.e., T = Ch"^, is taken when measurements are made in the norm, and a linear 
refinement path, i.e., t — Ch, when measurements are made in the norm. The 



results of an Cauchy error analysis are found in Tab. 5.3 and confirm second-order 
convergence in this case. The results of an Cauchy error analysis are found in 
Tab. |5.4]and confirm first-order convergence in this case. 



Our final test is a simulation of spinodal decomposition with different values of 



7. Specifically, we solve the scheme (5.1 )-(5.3) with = 0, i = 1, 2, 3, and with three 



values of 7; namely, 7 = 0, which yields the familiar Cahn-Hilliard model; 7 — 0.01; 
and 7 = 0.04. Furthermore, we take e = 0.01, h = r = 1 x 10"^, and T = 0.1. We 



use the same randomized initial data for the three simulations represented in Fig. 5.1 
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Y = 0.04 



Fig. 5.1. Spinodal decomposition for three values 0/7. The domain is Q = (0,1) X (0,1) and 
e = 0.01. The initial data are exactly the same for the three sim ulatio ns. The time step size is 
T = 1.0 X 10~^, and h = ^/256. We use a uniform mesh, as in Fig. The corresponding energy 

plots are shown in Fig. |5.2[ The average value of <p for all three simulations is approximately —0.1. 
For 7 = 0.01, the mass variation over the simulated time is only 1 X 10~^^. The max and min 
values of ip are very near the values +1 and —1, respectively. 



where the average value of ip is approximately —0.1. As expected, the mixture phase 
separates into domains wherein ip ^ —1 and ip ^ +1. Afterwards the system coarsens, 
as larger phase regions grow at the expense of smaller ones. The energy for the three 
simulations is displayed in Fig. |5.2| A general trend emerges, where, at least at 
early times, the energy decreases faster and the coarsening process is appears to be 
accelerated as the excess surface tension 7 increases. This behavior is expected and 
was observed in similar finite difference calculations undertaken in [25] , 

Note that we have proved that (at the theoretical level) the energy is non- 
increasing at each time step. This is observed in our computations. In addition 
to this, the mass, i.e, iph dx, at the theoretical level is expected to be unchanging 
from one time step to the next. On the practical level, we observe very little mass 
variation. For example, for the 7 = 0.01 case depicted in Fig. |5.1[ where initially 
'Poh dx w —0.1, we observe mass variation of only 1 x 10^^^ over the whole of the 
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0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

Time 



Fig. 5.2. Energy plots for the spinodal decomposition simulations depicted in Fig. \5.1\ The 
parameters for the simulations are given in the caption of Fig. \5.1\ The energy is observed to 
decrease at each time step. The general trend, at least at early times, is that the energy decreases 
faster with increasing values of the excess surface tension 7. 



simulation time. Note that our multigrid iteration stopping tolerance is of the same 



order, namely, tol — 1 x 10 in (A. 35 1. 
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Appendix A. Nonlinear Multigrid Solver. In this appendix we give the full 
details of the nonlinear multigrid solver that is used to march the scheme in time. 
Suppose C is polygonal, and assume that Te, i ~ 0,1, . . . , L, is a hierarchy of 



nested triangulations of fl as suggested in Fig. 5.1 In particular, Te is obtained by 
subdividing the triangles of Ti-i into 4 congruent sub-triangles. Note that hi-i = 2hg, 
i = 1, . . . , i, and that {Ti} is a quasi-uniform family. For simplicity, we shall use Pi 
finite element spaces and use the same space for the pressure as is used for the other 
variables. We define 

Vi^{ve C°(n) I v\k e Pi{K) yK e Ti} , 

for ^ = 0, . . . , L and observe the nested space chain Vq C Vi C V3 C • • • C Vl- Because 
of this nestedness, there is a natural injection operation Ig-i^i : Vg^i ^ Vg defined by 
Ii-i.i{v) — V, for all V G Vf-i, I = 1, . . . , L. Now, let Bi — {w£.i(x)}^-^ be the nodal 
basis for Vf, £ = 0, 1, . . . ,L. In other words, ug j (xii) = Sij, where are the 

nodes of Tg. We have level- wise representations of the unknowns of the form 

(A.l) ip({x) = '^ipe^iue,t{x) ■i==^ (p^ = {(pe^i,(pe^2, ■ ■ ■ ,'Pe,Ne)'^ , 

1=1 

and similarly for fJ-e(j^) and P£(x). Define the prolongation matrix via P^_i,f := \e-i,e, 
where is the Ni x A^^_i matrix representation of the injection operator Ii-i^i 

with respect to the bases and B^. There are two restriction operations — i.e., 

operations transferring information from the finer space to the coarser space V^-i 
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Fig. A.l. A hierarchical triangulation, Ti, = 0, 1, . . . , L, of a square domain O. Here L = 4, 
though in typical calculations we may use L = 8 or 9. 



— that we shall use. The first is called the canonical restriction and, in matrix form, 
is the Ni^i X Ng matrix defined via R^,^_i :— \J_i g The second is defined via 



AT;- 



(A.2) 



where the points X£_i_i are the nodes of the mesh Te-i- Note that {x^_i.i}^-^ ^ C 
{^i,i}f=i by construction. By Rg/-i we denote the matrix representation of Rij^i 
with respect to the bases Bg and Bg-i. 

In the present framework, our nonlinear finite element scheme is defined on the 
finest level, £ — L, as follows: find the triple {pL, HLj^l) & Vl x Vl x Vl such that 

(A.3) (VpL + 7<-1Vml, Vql) =0 Vgi e 14 , 

(A.4) +Te(VML,Vi^L) 

(A.5) [p.L.i>L)-e{^^^L,^i,L)-\({^LfvL,^l)^-\{vT\i'L) "^^L ^Vl. 

where t/'™ G is given. We have dropped the superscript m (the time step index) 
on the unknowns for simplicity. Theorem |3 . 4| guarantees that this problem always has 



a unique solution. The nonlinear system (A.3)-(A.5) may be written as 



(A.6) 
(A.7) 

(A.8) 



AlPl + tClMl = 0, 

e e 



where A^, B^, C^, M^, and Ql ((Pi) are x A^^ matrices whose components are 

(A.9) [Al],,, := (V^iL.^Vi^L,.), ■■= ([^T^f ^ u^.,,^ u^.^ , 

(A.IO) [Cl],,^. := (^r'V^iL,,, Vi^L,.) , [Mi],_^. := , 

(A.ll) [QLiVL)kj ■■= [i^LfuL,„UL,^). 



We solve (A.6|-(A.8) using a nonlinear multigrid method Ch. 5, §6]. This 
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requires that we split the equations into source (s) and operator (N) terms: 

(A.12) 0, (^J A^p^ + 7^/^^ , 

(A.13) sf:=Mi<^'»-\ (0J :=MiC^^+T(eAi + 7Bi)/x^+TCiPi, 

(A.14) := ^Micp^'^i, {ct,^)■.= e^LVL + -QL{^L)^L-^L^JiL. 

where := [p^, /i.^, y'/^] is the Nl x 3 array of unknowns. We must also define a 
"consistent" version of the nonlinear operator on all of the coarser levels. There are 
a number of ways to proceed in this task |6j; we choose the following path. Suppose 
that I £ {0, 1, . . . , L — 1} is given. We restrict the known solution from the previous 
time step to the coarser levels via 

/A ^ r\ I I D ^m— 1 / m— 1 m— 1 m— 1\ 

(A.15) iff, := 11 Rjj-iVl = ["Pe,! ^Vl2 ^■■■Vi,n,) > 

j=i+i 

Nt 
i=l 

Now, given any -0^ g Vi with the representation 

(A. 16) = {ilje,i,i^e,2, ■ ■ ■ ,i'i,Nt)'^ -^=> ^^(x) = ^ ■0^,iU£,i(x), 

i=l 

we define 

(A.17) [A,]^.^. (V?/,j, Vu,,,) , [B^]^_^. := (((^7"')' V?i,j, V^i,,,) , 

(A.18) [Q],,^^. (v^^-iVu,,,, Vm^,,) , [M,]^_^. := (m,,,, u,,,) , 

(A.19) [Q, := ((t/;,)' t/^,,, u,,,) . 

Observe that 

(A.20) = R^+i,fA£+iP^,f+i, M(,^Ri+i^iMi+iPe,e+i, 

which is standard in the finite element setting [6l [7] and is the reason for the term 
"canonical" describing Rg^ij. On the other hand, 

(A. 21) Bi « Rf+i^£B^+iPf^£+i, Ce « Ri+i,e^e+iPi,e+i- 

(Note that we could have recursively defined = R^+i^^B^+iPf^^-i-i, and similarly for 
Ci. But it turns out that this is an unnecessary complication from the point of view 
of the convergence of the algorithm.) Finally, we have 

(A.22) Nf)(C,) A,q, + 7Qiy,, 

(A.23) Nf ^ (C^) MixPf, + t {eAi + -/Be) + rQq^ , 

(A.24) {is) ekt^^ + ^Q^ (-0^) -0^ - M^z/^, 

where := [q^ , f ^, ■0^] is any given Ng x 3 array of unknowns. 
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We are now in a position to define the recursive nonlinear multigrid V-Cycle 
operator [HI Ch. 5, §6], which is the heart of our solver. In the following the superscript 
k is the V-Cycle loop index (not the time step index). Let </>^~^ := [p^"^, M^"^' V^"^] 
denote the current, level-£ multigrid iterate. For any Ni x S array of unknowns 



define Ni (|^) := 



N 



(1) 



these last two objects are Nix3 arrays by design. We define the action of the recursive 
nonlinear multigrid V-Cycle operator 



,(3) 



and 



„(1) „(2) J3) 



Note that 



(A.25) 

in the following 3 steps: 
1. Prc-smoothing: 
• Given (t)\~^ , 

(A.26) 



0^' = NMGM 



S£, A 



compute a smoothed level-£ approximation cf)f 



where 5 is a smoothing (or relaxation) operator, and A > is the number 
of smoothing sweeps. 
2. Coarse-grid correction: 

• Compute coarse-level initial iterate: 



(A.27) 



R 



• Compute the coarse-level right-hand side: 
(A.28) s^_i = Ri^i-i {si - N^(0^)) - 



(0£-i) 



• Compute an approximate solution of the following coarse grid equa- 
tion: 



(A.29) 



Note that this equation is uniquely solvable by Theorem |3.4| 
— If £ = 1 employ A smoothing steps: 



(A.30) 



(00. No, So) 



— If £ > 1 get an approximate solution to Eq. (A.29) using 4>(_i as 
initial guess: 

(A.31) = NMGM {£ - 1, N,_i,s,_i, A) . 

• Compute the coarse-grid correction: 
(A.32) 

• Compute the coarse-grid-corrected approximation at level k: 
(A. 33) (t>^ = Pc_i,f0£_i -f 4>f^. 

3. Post-smoothing: 
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Finally, compute 4>i by applying A smoothing steps: 



(A.34) 



cPf = S' 



When 



(A.35) 



37V, 



3 Nl 



(Pl 



< tol 



we stop iterating and set ^ 4>l — [PLjMliVl]; the fine-level solution. For 
smoothing, we use a nonlinear block Gaufi-Seidel method, like that discussed in |28j for 
a similar finite-difFerence nonlinear multigrid method. The exact details are omitted 
for brevity, but the principal idea is that the nodal values (pg^i, fi£^i,ipe,i) are always 
obtained simultaneously in the smoothing operation. We use A = 2 or 3 in the 
smoothing step. 
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